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ABSTRACT 

We compare the DARKexp differential energy distribution, N(e) oc exp(0o — e) — 1, 
obtained from statistical mechanical considerations, to the results of N-body simulations 
of dark matter halos. We first demonstrate that if DARKexp halos had anisotropic 
velocity distributions similar to those of N-body simulated halos, their density and 
energy distributions could not be distinguished from those of isotropic DARKexp halos. 
We next carry out the comparison in two ways, using (1) the actual energy distribution 
extracted from simulations, and (2) N-body fitting formula for the density distribution 
as well as N(E) computed from the density using the isotropic Eddington formula. 
Both the methods independently agree that DARKexp N(E) with (fro ~ 4 — 5 is an 
excellent match to N-body N(E). Our results suggest (but do not prove) that statistical 
mechanical principles of maximum entropy can be used to explain the equilibrated final 
product of N-body simulations. 

1. Introduction 

In this series of papers we revisit the use of statistical mechanics to address the structure of 
self-gravitating N-body systems, such as dark matter halos. In Hjorth Sz Williams (2010; Paper I) 
we showed that if it is correct to think of finite self-gravitating collisionless isotropic systems as the 
most probable configuration in the energy state space, then their energy distribution is given by a 
truncated exponential differential energy distribution, DARKexp; N(e) oc exp(<fio — e) — 1, where 
e and (fro are dimensionless energy and potential depth, respectively. The resulting structures have 
central density cusps. 

In Williams & Hjorth (2010; Paper II) we used the Extended Secondary Infall Model (ESIM) to 
confirm that restricted dynamical evolution can drive a system to the DARKexp state. ESIM largely 
fulfills the stipulations imposed by DARKexp; it is spherically symmetric, and its dynamics make 
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for efficient energy redistribution among particles, but does not allow the redistribution of angular 
momentum, which is not addressed in deriving DARKexp. ESIM generates a wide range of potential 
depths, as quantified by DARKexp's sole parameter, <j)Q. ESIM halos are not exactly isotropic, many 
show radial anisotropy, and some have tangential anisotropics. The anisotropy profiles vary between 
halos. These variations appear to be random; with no consistency, or 'universality' between halos. 

In contrast, N-body halos do show well defined anisotropy, which appears to be as universal as 

. The halo centers are isotropic but 
the outer parts have significant radial anisotropy. It is not surprising that N-body and ESIM halos 
have different anisotropy characteristics. ESIM dynamics include radial forces and accelerations, 
but no tangential accelerations (only velocities). Simulations, on the other hand, operate in full 
three dimensions, so tangential forces can redistribute particles' angular momenta and establish a 
universal anisotropy profile. Therefore it is possible that the presence of the well defined, universal 
anisotropy profile of N-body halos indicates an important element of the equilibrium structure of 
simulated halos, not accounted for in DARKexp models. 

A full theory for the distribution of particles in energy-angular momentum space using the 
maximum entropy principle is yet to be developed. In the meantime, we ask if we can limit 
ourselves to comparing energy distributions only, without regard to non-zero anisotropy. Since 
whatever produces the anisotropy profile in simulated halos may also effect the energy distribution, 
it is not immediately obvious that DARKexp N(E) can be compared to that of N-body halos. To 
proceed one must first determine what influence the anisotropy has on DARKexp N(E) and density 
profile. If the anisotropy of N-body halos is such that it has no effect on N(E) and p(r), then one can 
compare the DARKexp and N-body N(E) without worrying that anisotropy — and by implication 
whatever physics is present in simulations but was neglected in deriving DARKexp — will invalidate 
the comparison. 

In Section [2] we show that if anisotropy is similar to what the N-body simulations produce, 
then N(E) and p(r) are unchanged compared to their isotropic counterparts. Therefore we argue 
that a comparison between DARKexp and N-body N(E) can be justifiably made, and carry it out 
in two ways, as described in Sections |3 . 1 1 and |3 . 2| 

2. N(E, L) distributions of DARKexp models 

2.1. Generating self-consistent sets of N(E,L), p(r) and /3(r) 

A geometrically spherically symmetric non-rotating system is fully described by the distribution 
of its particles in the two dimensional plane of energy, E and angular momentum, L; N(E,L). 
The range of energies is bounded on one end by the value of the deepest central potential, and the 
escape energy, E = on the other end. For every E there is a maximum L which corresponds to 
the circular orbit at that energy; L max (E) is an upper envelope in the (E, L) plane. The lower 



the density profile (Hansen Moore]|2006 Hansen et al.||2010 



- 3 - 



boundary is L = 0, implying a radial orbit. An example of an L vs. E plot and the L max {E) 



envelope can be found in Wojtak et al. (2008), for the particular case of numerically generated 
universal halos. 

Even though the N(E, L) distribution contains all the information about a system, there is no 
simple relation (that we know of) between it and the corresponding density and anisotropy profiles. 

For the purposes of this work we need to check if a set of N(E,L), p(r), and anisotropy 



distributions is self -con sistent (Sections |2.3| and 2.4), or, generate self-consistent sets of these dis- 

2 1 \ 

tributions (Section 2.5). Anisotropy is defined as usual, /3(r) = 1 — ^j| r j i i n spherical systems the 
two tangential velocity dispersions are equal, ag = a^. We use the following procedure. Given 
input distributions N(E, L) and p{r) we numerically generate a halo by drawing particles from 
the N(E, L) distribution and adding up, or superimposing their orbits. The orbit superposition is 
done as follows. The system's potential <3?(r) is calculated from the input density profile. For each 
particle picked randomly from N(E, L) we find its radial velocity, v ra d(r), and apo- and peri-centers 
using the energy equation, E = $(r) + \ vf an (r) +v^, ad {r) , where the angular momentum is related 
to the radially dependent tangential velocity through vtanij") = L/r. The density that a particular 
particle contributes at a given radius is proportional to the amount of time it spends there, i.e., 
^ v rad d T '/ Jr ap ° V rad Thus a halo is built up from its individual particles. If the density profile 
generated in this fashion matches the input p(r), then the set of p{r) and /3(r) (computed afterward 
from the radial and tangential velocity dispersions) represent a self-consistent solution for the input 
N(E,L). 

In this work, the only input energy distributions and density profiles we use are those of the 
isotropic DARKexp models. In Sections |2.3| and |2.4| we use the above procedure to determine if 
specific L-distributions are consistent with DARKexp. In Section [2.5| we seek, through trial-and- 
error, an L-distribution that produces a specific anisotropy profile. So, if for a given L-distribution 
the input p(r) and the output p(r) from orbit superposition do not match then the L-distribution is 
rejected. (Note that accept ance/rejection is not done on individual particle basis; only whole halos 
are accepted or rejected.) The procedure is repeated until an acceptable L-distribution is found. 

This method has some features in common with the Schwarzschild orbit superposition method 
( Schwarzschild|1979 ) , where a self-consistent equilibrium halo is built up from its constituent orbits. 



There is no dynamical evolution. 



2.2. Anisotropy and the distribution of L 

To generate halos, one needs an input N(E,L). While DARKexp gives us the distribution in 
E, its says nothing about how orbits should be distributed in L to produce a given /3(r). To guide 
us, we start with some well known results connecting the energy distribution and the distribution 
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function, /. In general, 



N(E, L) dE dL = 8vr 2 L f(E, L) T r (E, L) dE dL, 



(1) 



where T r is the radial period of orbits (for example, see Appendix A of Wojtak et al. ( 2008[ >). For 
isotropic systems, / = f(E). Henon ( |1959 ) proved that only three types of potentials have T r that 
depends on E alone: the point mass, the uniform density sphere, and the isochrone potential (see 
also Efthymiopoulos et al. ( |2007 )). It is then more instructive to write eq. [T]as 



N(E, L 2 ) dE dL 2 = 4vr 2 f(E) T r (E) dE dL 2 



(2) 



The right hand side — and hence the left hand side — has no explicit L-dependence, therefore the 
number of orbits at a given E, between two L values is proportional to dL 2 , which means that for 
these three potentials isotropy implies a uniform distribution in I? . 

With this result in mind, one can rewrite eq. [T]to describe systems where orbits are distributed 
uniformly in L s , 

87T 2 



N(E, L ) dE dL° 



L 2 - s f T r dE dL". 



(3) 



If the distribution in a given parameter is uniform, that means the function has no explicit depen- 
dence on it. In the case of eq. [5J the expression for N(E, L s ) has no explicit L s , or L dependence. 
Further assuming T r = T r (E), implies that L 2_<5 / = h(E) is a function of E only, and hence / 
is separable; / = h(E)L s ~ 2 . These are systems of constant anisotropy, (5 = 1 — 5/2 (Henon 1973). 
A f3 = system is recovered for 5 = 2. Even though T r is not strictly a function of E only for a 
general system, it appears to be a very good approximation for many potentials. Hence a uniform 
distribution in L s can be used to approximate constant f3 systems of arbitrary density profiles. 



2.3. Isotropic halos 

In this section we generate N(E, L) distributions for isotropic DARKexp models. Although 



we know from Henon (1959) that a strictly uniform distribution in L 5 with 5 = 2 cannot produce 
isotropic systems for DARKexp, we nevertheless try a uniform distribution of particles in L 2 . This 
is shown as black points in the left panel of Figure [TJ for DARKexp (j)Q = 4. The magenta curve 
represents the upper boundary, L max (E). We apply the orbit superposition procedure described 
above, with an isotropic DARKexp N(E) and p(r) as the input, shown as dashed black curves in 
the top panel of Figure [2} The resulting density profiles are shown as the red, magenta and blue 
solid lines for 4>q = 2, 4, 8 DARKexp models, respectively. These agree very well with the input 
p(r), making this a self-consistent set of isotropic DARKexp halos. 

We conclude that isotropic DARKexp halos (with DARKexp form for the N(E)) are very well 
described by a uniform orbit distribution in L s , where 5 = 2. The deviations from this, which we 
know have to be present, appear to be smaller than the precision required in the present paper. 
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2.4. Anisotropic halos; /3 m 0.5 



For comparison, in the right panel of Figure [T] we show the case where, for any given E, the 
distribution is uniform in L, i.e. 5 = 1, so by the arguments of Section 2.2 these systems should have 
(nearly) constant (3 = 0.5. The anisotropy profiles of the orbit superposition halos are presented in 
the bottom panel of Figure [3] as the red, magenta and blue solid curves for <f>o = 2,4,8 DARKexp 
models, respectively. They do, in fact, have /3 « 0.5. 

The DARKexp input density profiles are the as black dashed lines in the top panel, while the 
orbit superposition profiles are represented by red, magenta and blue lines. These do not match 
the input p(r), and therefore DARKexp N(E), which we used as the input energy distribution, is 
not consistent with constant anisotropy as large as /3 = 0.5. 

A comment is in order. Inspection of Figure [3] shows that the systems obtained by orbit 



superposition obey the anisotropy-density inequality introduced by An Sz Evans (2006) for the 



central halo regions, and extended to apply to all radii by Ciotti & Morganti (2010). It states that 
7 > 2/3 (recall that p(r) oc r~ 7 ( r )). Applied to our systems it means that if /3 = 0.5, the density 
slope cannot be shallower than 7 = 1 anywhere in the system. This slope is represented by a 
dotted line segment in the upper panel of Figure [3j Our three density profiles are steeper than this 
at all radii. As a consequence of the anisotropy-density inequality, the oscillations in the density 
slope (Paper II) are erased, and the asymptotic slope is attained at much larger radii than in the 
isotropic models. 



2.5. Anisotropic halos; /3(r) « /3 N -bod y sim(r) 

Here we repeat the orbit superposition procedure but aim to generate anisotropy profiles 
similar to those of N-body simulated halos. Since these are not constant, a uniform distribution in 
L s will not work. Through trial and error we came up with a Gaussian distribution in L 2 . The full 
distribution in the E vs. I? plane is, 

N(e, L 2 ) de dL 2 cx N(e) x exp {-([I 2 / L 2 max (e)} - + a 2 ]f / (2a^ de dL 2 , (4) 

where N(e) oc exp(</>o — e) — 1 is the DARKexp form. As in Papers I and II we use dimensionless 
energy and potential, e = f3xE, and 4>q = (3tQo, where $0 is the system's potential depth, and 
Pt is its (negative) inverse temperature. The constants (01,02,03) are (0.25, 1.00,0.20) for $0 = 2, 
(0.25,0.75,0.15) for O = 4, and (0.25,0.50,0.10) for O = 8. 

The density and anisotropy profiles are shown as red, magenta and blue solid lines in the two 
panels of Figure |4j The anisotropy profiles generated through orbit superposition (red, magenta 
and blue curves in the bottom panel) are quite close in shape to the universal anisotropy profile. 



The latter is plotted (dotted line) using the Einasto profile 7 = 2(r/r__2)° with a = 0.17 (Navarro 



et al. 2004) combined with the anisotropy-density slope relation, 7 = —0.8 — 5/3 (Hansen et al 



- 6 - 



2010). The thick portion of the dotted line represents the region where the 7-/3 relation can be 



trusted (Navarro et al. 2010) 



Each of the three orbit superposition density profiles (red, magenta, blue) matches its cor- 
responding isotropic DARKexp profile (dashed black) very closely, which means that DARKexp 
N(E) remains a good description for systems with anisotropies as large as those seen in simula- 
tions of universal halos. In other words, DARKexp N(E) and density profiles originally derived for 
isotropic systems, are also consistent with systems whose f3(r) profiles are isotropic at the center 
and radially anisotropic closer to virial radii. We conclude that we can compare DARKexp N(E) 
to that of N-body simulated halos, and safely ignore the anisotropy of the latter. 



3. Comparing DARKexp N(E) to that of simulated halos 

Having demonstrated that universal anisotropy is too small to have an effect on N(E) we now 
compare the DARKexp energy distribution to that of simulated universal halos, using two methods, 
described below. 



3.1. Using N(E) of simulated halos 

Here we determine the differential energy distribution of dark matter particles residing inside 
the virial sphere of simulated halos. For this purpose we use a sample of 36 cluster-size relaxed halos 
extracted from a z = snapshot of a iV-body simulation of a standard ACDM cosmological model 

. This halo sample has 

already been used for calculating the distribution function and testing its phenomenological model 

. Each halo contains from 5 x 10 5 to 5 x 10 6 
particles inside its virial sphere defined in terms of the mean overdensity, {p)/p c ~ 100, where p c is 
the present critical density. 

We calculate N(E) for each halo independently by counting particles in energy bins. Then we 
combine all profiles into one and evaluate the median profile and the dispersion within the halo 
sample (blue curve and light blue area in Figure [5]). To combine the halos, we scale the binding 
energies by V 2 = GM{< r_2)/ r -2 and the particle numbers by N{< r_2), where r_2 is the radius 
where 7(r_2) = 2. These scalings preserve scale-free similarities of the phase-space properties of 
halos including the differential energy distribution itself. They also diminish the differences between 
N(E) profiles of individual halos with varying virial masses and scale radii. The scale radii, r_2 ; 
were obtained by fitting the NFW profile to the density profiles of simulated halos. We emphasize 
that the results are independent of the assumption of a fitting formula and remain the same if the 
scale radii are calculated by finding a maximum of p(r)r 2 , which is a fitting-independent method 
of measuring r_2- 



(for details of the simulation and the halo catalogue see |Wojtak et al][2008l 



with radially changing anisotropy (Wojtak et al.|| 2008 
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Although our analysis is restricted to the virialized parts of halos, we expect to encounter 
particle populations which are not fully equilibrated, e.g., particles which have only recently entered 
the virial sphere. In order to separate the energy range of these particles we indicate in Figure [5] the 
value of the gravitational potential at the virial radius r v of all halos (see the black line and light 
gray rectangle showing the median value and the la scatter within the halo sample). All relative 
energies, &o — E, below this value are covered by particles moving along the orbits confined inside 
the virial sphere so that this energy range likely represents a fully equilibrated particle population. 
On the other hand, the upper energy range is partially populated by orbits extending beyond the 
virial sphere so that it likely includes a particle population which has not settled to equilibrium. It 
is worth noting that at the energy separating these two energy regimes, i.e., E ~ <3?(r v ), N{E) starts 
deviating from an exponential growth predicted by the DARKexp model. At less bound energies, 
E > 3>(r v ), the presence of unequilibrated particle population manifests itself in an exponential 
cut-off of N(E). 

We fitted the median profile of N(E) for A-body halos with the DARKexp model given by 
N/Nq = exp[— (3t(E — $o)] — 1) parametrized by the normalization Nq and the inverse temperature 
(3t- The fitting was restricted to the energy range associated with the fully equilibrated particle 
population, i.e., <3?o — E < 3Vg (see Figure^. It yields (3t = 0.91 V^, which leads to the dimension- 
less potential depth of the halos <j)Q = 4.6. This best fit <fio is consistent with the most favourable 
range of values resulting from comparing DARKexp model with the best analytical approximations 



of the universal density profile {4>q ~ 4 — 5; Section 3.2). Given this agreement in <fto, and the results 



in Figure [5] we conclude that the DARKexp fit (dashed magenta curve) is a very good match to 
the N(E) profile from simulations. 



3.2. Using p(r) fitted to simulated halos 



Several fitting formulae have been suggested to describe the density profiles of N-body halos. 



The original double power-law NFW form was introduced by Navarro et al. (1997) as the universal 



shape for cold dark matter halos. It was later pointed out that simulated halos are even better 



fit with Einasto profiles (Navarro et al. 2004; Merritt et al. 2006). The S&M profile (Stadel et al 



2009) was recently proposed as a very good fit to the inner most regions, r < r_2- 



As in Paper II, one can compare the density profiles of numerically generated halos to those of 
isotropic DARKexp models. In Figure [6] we plot the density slope 7 vs. the log of the radius. The 
original NFW fitting formula (short-dash black curve) is fairly similar to DARKexp (three solid 
back curves). The Einasto profile (short- long-dash red straight line) tracks DARKexp considerably 
better than NFW, and does so for more than three decades in radius. While the Einasto is a 
straight line in the figure, DARKexp is 'concave', i.e., its density slope changes less rapidly with 
increasing radius. The S&M profile (long-dash blue curve) also has some concavity, and in this 
sense, it curves to mimic the DARKexp in the inner halo region. It is interesting to note that as 
the successive fitting formulae become better descriptors of the N-body halos (NFW — >• Einasto — > 
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S&M) they also get closer to the shape of the DARKexp density profile. 

Since DARKexp (<j>Q ~ 4 — 5) and N-body density profiles are very similar, and we already 
showed in Section |2.5| that anisotropy can be ignored for our limited purposes, we can assume 
isotropy and obtain N(E) from the density fitting functions of N-body simulations. From p(r) we 
first calculate f(E), using the isotropic Eddington formula, and then iV(JS]Q 

The results are shown in Figure [7] for five Einasto profiles (solid blue curves), parametrized by 
the a index of Navarro et al. ( |2004[ ); a between 0.07 and 0.27, in steps of 0.05. In calculating these 



we used 7 radial decades of p(r) between 10 5 r_2 and 100 r_2, which contain most of the systems' 
mass. We then fitted Einasto N(E) with the DARKexp form, N/N = exp[-f3 T (E - $ )] - 1. The 
fitting yields Nq and fix- As before, the product Pt^o is the dimensionless potential depth <po, the 
only shape parameter in DARKexp. 

In Figure [7] the DARKexp fits are represented by dashed magenta curves, and the best fit (f>o 
are indicated in the plot. Among these five cases the best match occurs at a ~ 0.17, which is also 



the average of halos simulated by Navarro et al. (2004). For a < 0.1 the two distributions do not 



fit at all; and, perhaps not coincident ally, such a values are not encountered in simulations. 



4. Conclusions 



It is often said that N(E) is primarily determined by p(r), and visa versa (Binney & Tremaine 



1987; Merritt et al. 1989), and that j3{r) has very little effect on this. We have quantified this 
statement for the specific case of DARKexp N(E) and p(r). We find that if DARKexp halos 
had /3(r) shapes similar to those of N-body simulated halos, their density and energy distributions 
could not be distinguished from those of isotropic DARKexp halos. This means that even without a 
full theory of collisionless equilibrium N(E,L), which would presumably explain (5{r) of simulated 
halos, one can compare the DARKexp and N-body N(E), while ignoring non-zero anisotropy. We 
have carried out this comparison in two ways, using (1) the actual energy distribution extracted 
from simulations, and (2) N-body p(r) fitting functions, and N(E) computed from p(r) using the 
isotropic Eddington formula. Both of these methods agree that DARKexp N(E) with <po ~ 4 — 5 
is an excellent match to N-body N(E). This suggests that statistical mechanical principles of 
maximum entropy can be used to explain the equilibrated final product of N-body simulations. In 
the future we will extend the maximum entropy principle used to derive DARKexp to include L. 
This will hopefully explain f3(r) observed in simulated halos. 



1 One might be tempted to use the isotropic Jeans equation, thereby obtaining the velocity dispersion profile 
c(r), and then estimating kinetic energy as |<r 2 . Combining this with the potential gives E for particles at that 
radius. This procedure will give wrong results because the velocity distribution function (VDF) is implicitly assumed 
to be Maxwellian, which is wrong. The importance of VDF and its deviations from Maxwellian was discussed in 



Kazantzidis et al. (2004) and Hansen (20091. 
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Fig. 1. — The log of angular momentum vs. the energy of particles making up halos. The energy 
distribution is the DARKexp N(e), with <j>Q = 4 and particles shown are a realization of this 
distribution. The magenta line is the envelope corresponding to circular orbits, L max {E). Left 
panel: For a given E, the distribution in L 2 is uniform between and L 2 nax (E). This gives an 
approximately isotropic system (Section |2.2[ Figure[2]). Right panel: For a given E, the distribution 
in L is uniform between and L max (E). This gives constant j3(r) ~ 0.5 (Section |2.2[ Figure [3]). 
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Fig. 2. — Density and anisotropy profiles for DARKexp models (dashed black curves), and halos 
generated using the method described in Section [2~T] (solid color curves). To achieve isotropy for 
the red, magenta, and blue solid curve halos we used a uniform distribution in L 2 , illustrated in 
the left panel of Figure [T] 



- 13 - 



n 1 1 r 



"i 1 1 r 



1 1 1 r 



00. 



2 - 



I 



3 



o 



-2 - 



-4 
1 



0.5 







-0.5 



black dashed lines: 
DARKexp O =2,4,8 




1 
log (r/r_ 2 ) 



Fig. 3. — Density and anisotropy profiles for DARKexp models (dashed black curves), and halos 



generated using the method described in Section 2.1 (solid color curves). At every E, the distri- 
bution in L is uniform between and L max (E); this distribution is illustrated in the right panel 
of Figure [TJ and produces halos with /3 ~ 0.5. These systems obey the density-anisotropy relation, 
7 > 2/3 since the density slope is never shallower than 7 = 1, shown as the thin dotted line segment 
in the upper panel. 
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Fig. 4. — Density and anisotropy profiles for DARKexp models (dashed black curves), and halos 
generated using the orbit superposition method described in Section 2.1 These density profiles (red, 
magenta, and blue curves) agree very well with the isotropic DARKexp p(r). The distribution in 
the L vs. E plane (see Section 2.5) was chosen such that the anisotropy profile would look similar 
to that of N-body halos. The latter is plotted (dotted line) using the Einasto profile 7 = 2(r/r_2) Q 
with a = 0.17 (Navarro et al. 2004) combined with the anisotropy-slope relation, 7 = —0.8 — 5/3 
(Hansen et al. 2010). The thick portion of the dotted line represents the region where the 7-/3 
relation can be trusted. 
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Fig. 5. — Comparison between N(E) from iV-body simulation (blue curve) and the DARKexp fit 
with (J)q = 4.6 (dashed magenta curve). The blue curve and light blue area represent the median 
profile and the la scatter within the halo sample. The light gray area indicates the gravitational 
potential at the virial sphere (with the width corresponding to the la scatter within the halo sample 
and the black line showing the median value) . It defines an upper limit of the relative energy for a 
set of particle orbits which are fully confined inside the virial sphere. 
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Fig. 6. — Comparison of DARKexp density profile slopes (three solid black curves) with N-body 



fitting functions: Einasto (Navarro et al. (2004) a = 0.17, short-long-dash red straight line), NFW 



(short-dash black curve), and S&M (Stadel et al. (2009); long-dash blue curve). The last profile is 
meant to fit only the inner region of N-body halos, so the radial range outside of r_2 is not plotted. 
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Fig. 7. — Comparison between N(E) for DARKexp and N-body halos; see Section 3.2 A range of 
Einasto profiles (solid blue), with a = 0.07,0.12,0.17,0.22,0.27 were fit with DARKexp (dashed 
magenta), and the fitted (po values are shown in the plot. The best match occurs for the a value 
(~ 0.17) which happens to be the average of simulated halos (Navarro et al. 2004). For a < 0.1 
the two distributions do not fit at all, and these a values are not encountered in simulations. 



